Theorem (Dirichlet, 1842) Let $\alpha,Q\in\mathbb{R}$ such that $Q>1$. Then there exist integers $p,q$ such that $1\leq q<Q$ and $$\left|\alpha q-p\right|\leq {1\over Q}$$.
Corollary If $\alpha\in\mathbb{R}\setminus\mathbb{Q}$, then there are infinitely many pairs $(p,q)$ of coprime integers such that $$\left|\alpha - {p\over q}\right|<{1\over q^2}.$$
In [23]:
def diophantine_approx(alpha, Q):
for q in range(1, Q):
if frac(alpha*q) <= 1./Q or 1-1./Q <= frac(alpha*q):
p = round(alpha*q)
return (p,q)
In [4]:
%display latex # output in latex
In [26]:
diophantine_approx(pi, 10)
Out[26]:
In [27]:
n(22/7 - pi)
Out[27]:
In [28]:
diophantine_approx(pi, 100)
Out[28]:
In [29]:
diophantine_approx(pi, 200)
Out[29]:
Finding the next convergent takes time:
In [30]:
%time diophantine_approx(pi, 30000)
Out[30]:
In [31]:
%time diophantine_approx(pi, 35000)
In [2]:
cf_pi = continued_fraction(pi)
cf_pi
Out[2]:
In [5]:
latex(cf_pi[:25])
Out[5]:
In [35]:
L = cf_pi.convergents()
L[:20].list()
Out[35]:
In [36]:
numerical_approx(pi - 103993/33102)
Out[36]:
In [37]:
numerical_approx(1/33102^2)
Out[37]:
Theorem (Dirichlet, 1842) Let $\alpha_1,\dots,\alpha_d\in\mathbb{R}$ and $Q>1$ an integer. Then there exist integers $p_1,\dots,p_d,q$ such that $1\leq q<Q^d$ and $$|\alpha_i q-p_i|\leq {1\over Q}\quad\text{ for all } 1\leq i\leq d.$$
Corollary If at least one of the $\alpha_1,\dots,\alpha_d$ is irrational, then there exist infinitely many $(d+1)$-uplet $(p_1,\dots,p_d,q)$ such that $$\left|\alpha_i - {p_i\over q}\right| < {1\over q^{1+\frac{1}{d}}}\quad\text{ for all } 1\leq i\leq d.$$
Last month, I coded this (it will be available in the next release of slabbe package):
In [38]:
from slabbe.diophantine_approx import best_simultaneous_convergents
It uses Cython code and computations are done in parallel on all the available cpus:
In [39]:
from itertools import islice
In [41]:
it = best_simultaneous_convergents([e,pi])
%time [(p1,p2,q) for p1,p2,q in islice(it,18)]
Out[41]:
A Multidimensional Continued Fraction (MCF) algorithm is a function $$ \begin{array}{rrcl} F: & \Lambda & \to & \Lambda\\ & x & \mapsto & M(x)^{-1}\cdot x. \end{array} $$ where $\Lambda\subset\mathbb{R}^d$ is a cone and $M:\Lambda\to SL(d,\mathbb{N})$ is
Classical references are Schweiger 2000, Brentjes 1981.
The positive cone $\Lambda=\mathbb{R}^3_+$ is partionned as $\Lambda=\cup_{\pi\in\mathcal{S}_3}\Lambda_\pi$ where $$ \Lambda_\pi = \{(x_1,x_2,x_3)\in\Lambda\mid x_{\pi 1}< x_{\pi 2}< x_{\pi 3}\}. $$ Matrices are $M(\mathbf{x}) = M_\pi$ where $$ \begin{array}{lll} M_{123}={\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 1 & 1 \end{array}\right)} & M_{132}={\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 1 \\ 0 & 0 & 1 \end{array}\right)} & M_{213}={\left(\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 1 & 0 & 1 \end{array}\right)} \\ M_{231}={\left(\begin{array}{rrr} 1 & 0 & 1 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)} & M_{312}={\left(\begin{array}{rrr} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)} & M_{321}={\left(\begin{array}{rrr} 1 & 1 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}\right)} \\ \end{array} $$ For example, if $(x_1,x_2,x_3)\in\Lambda_{123}$, then $$ F (x_1,x_2,x_3) = (x_1, x_2, x_3-x_2). $$ In other words, it subtracts the second largest entry to the largest.
In [42]:
from slabbe.mult_cont_frac import Brun
F = Brun()
In [43]:
orbit = F.cone_orbit_list((10, 21, 37), 5)
[(x,y,z) for (x,y,z,_,_,_,_) in orbit]
Out[43]:
In [ ]:
orbit = F.cone_orbit_list((e, pi, 1), 5)
[(x,y,z) for (x,y,z,_,_,_,_) in orbit]
The algorithm $F$ defines a cocycle $M_n:\Lambda\to SL(d,\mathbb{Z})$ $$ M_0(\mathbf{x}) = I \quad \text{and} \quad M_n(\mathbf{x}) = M(\mathbf{x}) M(F\mathbf{x}) M(F^2\mathbf{x})\cdots M(F^{n-1}\mathbf{x}). $$ with the cocycle property $M_{m+n}(\mathbf{x})=M_m(\mathbf{x})M_n(F^m(\mathbf{x}))$.
In [44]:
[F.n_matrix((e,pi,1), i) for i in range(5,25)]
Out[44]:
In [45]:
from slabbe.diophantine_approx import best_simultaneous_convergents
from slabbe.diophantine_approx import mult_cont_frac_vs_dirichlet
from slabbe.mult_cont_frac import Brun,ARP,Reverse,Selmer,Cassaigne
In [46]:
v = [e, pi]
it = best_simultaneous_convergents(v)
dirichlet = [next(it) for _ in range(18)]
algos = [Brun(), ARP(), Reverse(), Selmer(),Cassaigne()]
mult_cont_frac_vs_dirichlet(v, dirichlet, algos)
Out[46]:
In [ ]:
To a MCF algorithm corresponds a dynamical system $(\Delta, f)$ defined by the projection of the map $F$ on a simplex $\Delta$: $$ f(\mathbf{x}) = \frac{F(\mathbf{x})}{\Vert F(\mathbf{x})\Vert} \qquad\text{ on }\qquad \Delta=\{\mathbf{x}\in\Lambda\mid\Vert\mathbf{x}\Vert=1\}. $$ Below we illustrate some examples coming from :
Now we illustrate how to reproduce experimentations available in http://arxiv.org/abs/1511.08399
using slabbe Sage Optional Package :
sage -pip install slabbe
In [52]:
from slabbe.mult_cont_frac import Brun
F = Brun()
cocycle = F.matrix_cocycle()
t = cocycle.tikz_n_cylinders(3, scale=3)
print repr(t)
In [50]:
t.pdf()
Out[50]:
In [53]:
%display default
%matplotlib inline
In [54]:
from slabbe.mult_cont_frac import Brun
F = Brun()
_ = F.invariant_measure_contour_plot(1000000, 40, norm='1')
_ = F.invariant_measure_wireframe_plot(1000000, 40, norm='1')
In [55]:
F.lyapunov_exponents(n_iterations=10^8)
Out[55]:
In [56]:
from slabbe.lyapunov import lyapunov_comparison_table
In [57]:
from slabbe.mult_cont_frac import Brun,ARP,Reverse,Selmer,Cassaigne
from slabbe.lyapunov import lyapunov_comparison_table
algos = [Brun(), ARP(), Reverse(), Selmer(),Cassaigne()]
lyapunov_comparison_table(algos, n_iterations=10^5)
Out[57]:
In [58]:
V = [[1,0,1],[1,0,0],[1,1,0],[0,0,-1],[0,1,0],[-1,0,0],[0,1,1],[0,0,1],[0,-1,0]]
P = Polyhedron(vertices=V).polar()
P
Out[58]:
In [60]:
s = P.projection().tikz([674,108,-731],112)
In [62]:
type(s)
Out[62]:
In [ ]:
view(s) # slow, sometimes broken...
In [63]:
from slabbe import TikzPicture
t = TikzPicture(s)
t
Out[63]:
In [67]:
from IPython.display import Image
In [72]:
Image(t.png(view=False))
Out[72]:
In [73]:
sage: S = FiniteSetMaps(5)
sage: I = S((0,1,2,3,4))
sage: a = S((0,1,3,0,0))
sage: b = S((0,2,4,1,0))
sage: roots = [I]
sage: succ = lambda v:[v*a,v*b,a*v,b*v]
sage: R = RecursivelyEnumeratedSet(roots, succ)
sage: G = R.to_digraph()
sage: G
Out[73]:
In [74]:
t = TikzPicture.from_graph(G)
In [75]:
Image(t.png(view=False))
Out[75]:
In [76]:
TikzPicture.from_graph?
In [78]:
TikzPicture?
In [ ]: